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Abstract 

Monte Carlo simulations are presented for a coarse-grained model of real quadrupolar fluids. 
Molecules are represented by particles interacting with Lennard-Jones forces plus the thermally 
averaged quadrupole-quadrupole interaction. The properties discussed include the vapor-liquid 
coexistence curve, the vapor pressure along coexistence, and the surface tension. The full isotherms 
are also accessible over a wide range of temperatures and densities. It is shown that the critical 
parameters (critical temperature, density, and pressure) depend almost linearly on a quadrupolar 
parameter q = Q*^/T*, Q* is the reduced quadrupole moment of the molecule and T* the reduced 
temperature. 

The model can be applied to a variety of small quadrupolar molecules. We focus on carbon 
dioxide as a test case, but consider nitrogen and benzene, too. Experimental critical temperature, 
density and quadrupolar moment are sufficient to fix the parameters of the model. The result- 
ing agreement with experiments is excellent and marks a significant improvement over approaches 
which neglect quadrupolar effects. The same coarse-grained model was also applied in the frame- 
work of Perturbation Theory (PT) in the Mean Spherical Approximation (MSA). As expected, the 
latter deviates from the Monte Carlo results in the critical region, but is reasonably accurate at 
lower temperatures. 
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I. INTRODUCTION 



Solvents play an essential role in the design and processing of many molecular materials 
(e.g., oligomers, poljmaers, etc). In comparison to a melt, the molecular mobility of dissolved 
substances increases considerably in solution. Not only the flow properties can be controlled 
easily in solution, but also the phase behavior (and hence the morphology) of the dissolved 
materials, e.g., by changing the thermodynamic state conditions like temperature, pressure, 
and concentration. 

A particularly important solvent is supercritical carbon dioxide, because the material is 
inexpensive, nonpoisonous, not reactive, and thermally stable. Hence, its application as a 
solvent is widespread.-»2i'^ However, the phase behavior of polymer- solvent systems or other 
binary fluid mixtures is rather complex in general. When the thermodynamic control param- 
eters temperature T, pressure p and solute molar fraction x are varied, various liquid-vapor 
and fiuid-fiuid phase equilibria occur, and many different types of (rather complicated) phase 
diagrams can be observed^"^. Even for simple binary fluid mixtures, e.g., carbon dioxide plus 
short alkanes such as hexadecane, the phase diagram is only known rather incompletely from 
experiment.-*^ These uncertainties also hamper the judgment of the accuracy of the theoret- 
ical modeling of such systems.-'^'i^i'ii In fact, due to the large control parameter-space that 
needs to be explored, comprehensive experimental work would be very cumbersome, and a 
modeling approach seems to be the method of choice. However, the large number of states 
(T, p, x) that need to be simulated and the complexity of the systems renders a fully chemi- 
cally realistic all-atom simulation practically impossible. Thus, the construction of a suitable 
coarse-grained model for such systems containing polymers (or oligomers, respectively) is 
very desirable. While there is a rich literature on the construction of coarse-grained models 
for (flexible) polymers,— ii^'i^>i^ii^>i^ii^>^ comparatively little attention has been paid to the 
question on how a coarse-grained solvent molecule such as CO2 should be described. Iwai 
et al.— and Virnau et al.-'^'^ simply used particles interacting with simple Lennard- Jones 
forces among themselves and with the beads of the bead-spring chain that represents effective 
subunits of the polymer. While particles interacting with Lennard- Jones potentials describe 
noble gases such as liquid argon or neon rather well, it is clear that a "Lennard- Jonesium" is 
a somewhat unsatisfactory description of a carbon dioxide molecule. While considerable at- 
tention has been paid to atomistic models of CQ^2^^^^^^^i2^i2L2^i^^^^iM.^i^ -we are 
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not aware of comprehensive systematic studies of coarse-grained models for this molecular 
fluid. With respect to the atomistic models for CO2, we note that there is no consensus in 
the literature on a unique form of the interaction potential and its parameters. Starting from 
the Murthy-Singer-McDonald (MSM) model,— several potentials have been proposed (for a 
recent comparison see Ref. i2S ). In Ref. |27|, two variants of the Elementary Physical Model 
(EPM) force field were suggested, that yielded critical temperatures of Tc = 313.4 ± 0.7 K 
and 312.8 ± 3.0 K, respectively, while the experimental value is Tc = 304.2 K.— In view of 
this 3% discrepancy between the atomistic models and the experiment, it was suggested^I to 
use the experimental critical temperature and rescale the energy parameters of the model to 
reproduce the correct value of the critical temperature (EPM2). In fact, Virnau et al.-*^*^, 
using a simple "Lennard-Jonesium" to model CO2, fixed the Lennard- Jones parameters to 
match both the critical temperature Tc and the critical density pc- Since the atomistic 
models underestimate the critical density (yielding^^ 453.7 ± 4.3 kg/m^ or 449 ± 16 kg/m^ 
instead of the experimental value^^ pc = 468.0 kg/m^), they also require a corresponding 
rescaling of the interaction range parameters. Hence, EPM2 needs the same input from 
the experimental critical data as the coarse-grained model of Virnau et al.-'^'^. For the 
resulting model, the coexistence densities predicted for the liquid branch in the temperature 
region 230K < T < 280K deviate distinctly less from the experimental results^- than the 
corresponding results of the coarse-grained modeL-^i^ 

As indicated above, the main interest for obtaining an accurate coarse-grained model for 
CO2 is its potential application in multicomponent systems, e.g., polymer solutions in which 
CO2 acts as a solvent.-*^ For such systems also many attempts were undertaken to derive 
approximate analytical equations of states (e.g., Refs. SljsS) and it is, of course, also highly 
desirable to validate such equation of state theories by simulations. However, the coarse- 
grained model for CO2 of Virnau et al.*''°, when combined with a suitable coarse-grained 
model for the alkanes, required rather large deviations from the simple Lorentz-Berthelot 
mixing rules to account for the available experimental data.-"^ Most likely, the somewhat 
oversimplified CO2 model is responsible for most of these deviations. Approximating CO2 as 
a Lennard- Jones particle without considering its rather large quadrupolar moment (IQI = 4.3 
DA) is probably not sufficient - the unit D (Debye) equals 10~^^ in CGS units which are 
adopted throughout the manuscript. 

In the present paper we explore a slightly more involved coarse-grained model for 002.—*^ 
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The molecule is still described as a Lennard- Jones particle, but we also include the exper- 
imentally known quadrupole moment as an input parameter, together with critical tem- 
perature and critical density. A precondition for the usefulness of coarse-grained models 
is that simulation codes execute very fast. The angular-dependent quadrupole-quadrupole 
interaction requires significant computational resources which would be a serious drawback 
to such a model. However, compared to the Lennard- Jones forces, the quadrupolar interac- 
tion is still a rather weak perturbation. Therefore, we apply one further approximation:— 
the angular dependence is averaged over in a second order thermodynamic perturbation 
calculation. Thus, an effective isotropic potential is obtained. Rather encouraging results 
using such an approximation have been reported in the literature. Miiller and Gelb^ es- 
timate coexistence curves from non-equilibrium molecular dynamics (NEMD) simulations 
of temperature quenches from the one-phase region into the two-phase region, where one 
then waits until the system has phase separated into the two coexisting phases.— In this 
manuscript we apply grand-canonical Monte Carlo methods,— i^i^ combined with a finite 
size scaling^'^'^ analysis. This allows us to locate precisely the critical point of the model. 
Note that a direct estimation of the critical point from the simulation is difficult if either 
Gibbs ensemble techniques'^^ '^^ or the temperature quench technique^i^ are applied. In 
these cases, one relies on a fit of the coexistence data to a suitable power law extrapolation. 
With the present techniques one can obtain the critical properties very accurately. This 
precision is required because the critical properties are used to gauge the Lennard- Jones 
parameters of the model. 

In Sec. II, we give a more detailed description of our model and simulation techniques. 
Sec. Ill describes our results for carbon dioxide and compares them to previous approaches. 
Sec. IV discusses the application of the model to other quadrupolar fiuids, namely nitrogen 
and benzene. Sec. V describes the application of first order perturbation theory in the 
mean spherical approximation (PT-MSA) to precisely the same model which was used in 
the simulation, thus allowing a meaningful comparison. Finally, Sec. VI concludes the 
discussion and gives an outlook on future work. 
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II. MODEL AND SIMULATION TECHNIQUE 



A. Choice of Model 



Our model system consists of neutral spherical particles which carry a quadrupolar mo- 
ment Q and interact with each other both via the Lennard- Jones potential 



U^'^ = Ae 



r,- 



and the quadrupole-quadrupole interaction^ 



a \6 



(1) 



(2) 



The angle- dependent part is given by: 



/^^ = 1 - 5 cos^ Oi-b cos^ Oj + 17 cos^ Oi cos^ % 
+2 sin^ 9i sin^ 9j cos^(0j — (pj) 
—16 sin 6'j cos 9i sin 6'j cos 6'j cos(0i — (pj) . (3) 

In Eqs. ([H E]), = |rj — -r^l is the distance between molecules at sites rj, rj, while (6'j, 0j) 
are the polar angles characterizing the mutual orientations of the (linear) molecules ( 6i are 
the angles between the axis joining the two molecules and the quadrupole vectors of the 
molecules; (pi are the rotational orientations of the quadrupole vectors relative to the joining 
axis). In Eq. ([T]), e and a set the scales of energy and distance for the Lennard- Jones (LJ) 
interaction, respectively. 

The angular-dependent part of the potential (Eqs. ([2]), ([3])) slows down the speed of 
the algorithm considerably. Therefore, following Ref. l39j, we average over the angles of the 
quadrupolar potential to create an effective isotropic representation. More precisely, one 
expands the Boltzmann factor exp{-~j3Ufj^), {(3 = (fc^T)"^), in a Taylor series to second 
order in /3. After taking averages over the angles, the following temperature-dependent 
isotropic potential is obtained: 



For the potentials of Eq. ([T]) to Eq. (jlj) one can employ the standard procedure^i^ of 
cutting and shifting to zero at a cutoff distance rij = Vc = 2\/2a typically applied to 
Lennard- Jones systems. The total potential then reads 



12 / \ 6 / \ 10 
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As 

, r > Vc- 

The reduced quadrupolar interaction parameter is defined as 



r < Tr. 

(5) 



and Tc are the values of the reduced quadrupole parameter and temperature at the critical 
point. 5* shifts the cut potential to zero at r^j = r^, so that U{rij) is continuous everywhere: 

S=^+^-^ (7) 
16384 5 256 ■ ^ ^ 

Note that Eq. ([6]) is given in CGS units. In SI units, there is an additional factor (47r£:o)^^- 
It is clear that U{rij) is explicitly temperature-dependent because q and S are temperature- 
dependent. Hence, special care needs to be exerted when temperature derivatives are taken. 
For instance, the fluctuation relation linking the specific heat to the fluctuations of the 
potential energy no longer holds for Eq. ([5]). We also note that Eq. ([5]) differs from the 
potential obtained when one cuts off Eqs. (|T|2|) at rij = Tc- Indeed, continuity of U^^ 
would require an orientation-dependent shift of the potential. It is also well-known^^ that 
the relation between the critical temperature of a fluid and the energy scale e of the LJ 
interaction depends rather strongly on the cutoff r^. Our choice of a rather small value 
for the cutoff is mainly motivated by the desire to have a very fast simulation algorithm, 
but larger cutoffs will lead to very similar results. As we will demonstrate later, differences 
in the phase diagram almost disappear when simulation data are rescaled to match the 
experimental critical point for different r^. A further motivation for this choice of the cutoff 
is that for g = our model reduces to that of Refs. 

Our strategy will be to compute the critical temperature Tc{qc) and the critical density 
pciQc) from the simulation, using the potential from Eq. ([5]). Following previous work-»^, 
e and a are determined by the condition that these critical parameters match precisely 
their experimental counterparts. In the following, T* and p* will refer to temperatures 



and densities (and other quantities that will be introduced with "*") expressed in units of 
e{qc), cr{qc) and MmoI, the molar mass of the fluid. We need to consider that the parameter 
Qc, Eq. depends itself on e and a, and not only on the (given) experimental value for 
Q. This difficulty is related to the quadrupolar interaction in Eq. which shifts the 
critical point in the (T, p) plane relative to its position for Q = 0. Even if one is only 
interested in a single choice of Q, a simulation of a single model system (i.e., one choice of 
e, a and Q or g, respectively) is never sufficient to deal with this problem. However, this 
puzzle can be solved by determining the critical lines T*{q^ and p*{qc) as a function of the 
(dimensionless) parameter qc- Fig. [1] shows the results of this calculation and demonstrates 
that both T*{qc)/T*{0) and pl{qc)/Pc{^) are very smooth functions of qc- These curves are 
almost linear, so recording a few (altogether 9) choices of nonzero qc was sufficient to obtain 
good accuracy. In the range < qc < 0.5, the critical temperature increases by almost 30% 
while the critical density increases by about 10%. 

Having determined T*{qc) and Pciqc), one can compute easily e{qc) and o"(gc) such that 
the model corresponds to a specific experimental system Tc^exp and Pc,exp- Eq- (E]) must hold 
together with 

e{qc) = ksTc,..,/T:{qc), a3(ge) = [^^^^^^1 • (8) 

Here Mmoi is the molar mass of the simple molecule and A^^ Avogadro's number. These 
equations are solved by a simple iteration procedure, using the following fit functions repre- 
senting the data of Fig. [H 



Tc{qc)/T*{0) = 1 + 0.46111 g, + 0.17571 (9a) 
P:fe)/p:(0) = l + 0.19298ge , (9b) 

where T*{qc = 0) = 0.99821 and p*{qc = 0) = 0.32276. Appendix [B] explains in detail how 
simulation parameters are derived from experimental data. We note that the limiting factor 
for the accuracy of our procedure is not at all the limited accuracy of Eqs. ( l9al [9bl) . but 
rather the uncertainty with which the physical quadrupole moment Q of the molecule, needed 
as an input to Eq. ([6]), is known. Considering CO2 as an example, we take Q = (—4.3 ±0.2) 
DA. However, since Q is raised to the fourth power in Eq. ([6]), the 5% uncertainty in Q 
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becomes a 30% uncertainty in the reduced simulation parameter q. For Q = —4.3 DA, we 
obtain 

qc = 0.387, e = 3.491 x lO^^V, a = 3.785 A (10) 

The uncertainty in Q would actually allow for a range 0.32 < q^ < 0.47 with corresponding 
changes of e and a. In view of these uncertainties, one could not hope for a perfect agreement 
between the simulation results (for other quantities rather than pc and Tc) and experiment, 
even if the form of the coarse-grained potential, Eqs. ([MZI), were perfectly accurate. 

Already at this point, we note that nothing in the model (Eqs. ([5])- ([7])) is specific to CO2. 
Hence, Fig. [T] (or Eqs. fl^ [9b|) . respectively) can be used for modeling other quadrupolar 
fluids, too. This fact will be taken up in Sec. IV and Appendix [Bl We also note that e and 
a are independent of the state of the system once they are fixed, g, however, is given by 
q = qc ■ Tc/T (according to Eq. IQ), which needs to be considered when coexistence curve 
and interfacial tension are calculated. 

B. Comments on the Simulation Technique 

In this section we comment briefly on the Monte Carlo simulation techniques which 
are required for the computation of Fig. [T] and other physical properties. As in previous 
work,-ii^ extensive simulations were undertaken in the fiVT ensemble, where the box volume 
V = L^, the chemical potential fi of the particles and the temperature are fixed. The particle 
number fluctuates, since the elementary Monte Carlo move consists of random insertions or 
deletions of particles. Thus, long wavelength fluctuations of the density are equilibrated 
easily. In contrast. Molecular Dynamics or canonical ensemble Monte Carlo methods that 
conserve the particle number in the system suffer from a slow equilibration of long wavelength 
density fluctuations ( "hydrodynamic slowing down"^'*). The temperature quench simulations 
encounter the additional difficulty that vapor-liquid interfaces extending throughout the 
simulation box are formed. Such interfaces are notoriously slowly relaxing and strongly 
fluctuating objects and thus avoided in Gibbs ensemble techniques.— 

For the sake of efficiency, histogram extrapolation techniques are used. In a typical MC 
run, the particle number n and the total energy E are recorded at regular intervals. The 
resulting distribution P^^xin^E) can then be extrapolated to neighboring values of /x' and 
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T' using the following expressio: 

P,'An,E) = ^PMn,E)exp [(f^ - f )^ " " ^)^] (H) 

with A/" being a normalization constant. Here, we have assumed that Qc remains constant. 
Extrapolations at constant Q would require an additional reweighting factor related to the 
temperature-dependence of the potential (Eq. [5]). Of course, Eq. ffTTl) is only accurate when 
P^^T{n,E) and P^'^T'{n, E) overlap strongly. Nevertheless, reweighting is very useful for fi 
near ficoex{T), where two-phase coexistence between vapor and liquid occurs. In this region, 
Pii,T{n) = J dEP^ T{n, E) has a two-peak structure: one peak occurs at pcoL ~ n/V^ the 
vapor density at coexistence, the other peak at pioL ~ n/V , the liquid density at coexistence. 
For yU = ficoex{T), the areas underneath both peaks are equal ("equal area rule"— i^), but 
unfortunately ficocxiT) is not known beforehand. However, if one has Pf^^T{n), for some n 
close enough to ficoexiT), one can try to reweight the data according to Eq. f|TT]) with no 
additional simulation effort. In this way, the coexistence curve can be located precisely. 
The corresponding pressure is computed from the virial equation. All these procedures have 
already been applied in previous work for = 0. For more details the reader is referred to 



Refs. 



IG 



Following a path along p = ficoexiT) in the {fi,T) plane and recording moments of the 
density distribution, we calculate 2^^*^ and 4**^ order cumulants 



f/2 = (M^)/(|M|)2 , U,= {M')/{MY , M^p-{p). (12) 

Reasonably accurate estimates for Tc can be obtained from the intersection point of 
either U2{T) or Ui{T) for different L. The justification of this simple recipe follows from 
the theory of finite size scaling.— i^i^i^i^ Fig. [2] shows that T^. can be determined with a 
relative accuracy of about 3/10'^ with moderate computational effort. The lack of perfect 
intersections in the size range 9(t < L < 13.5(7 indicates that the asymptotic region of 
finite size scaling has not been reached yet, and corrections to finite size scaling are still 
present. However, the estimate ksTc/e = 1.152 ± 0.003 is clearly accurate enough for our 
present purposes. Note that the simple analysis presented in Fig. [2] ignores "field mixing" - 
effects^ between density and energy per particle. Of course, for a high precision study of 
critical exponents and critical amplitudes, more sophisticated finite size scaling methods are 
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available,— but this is beyond the scope of the present investigation. 

For temperatures distinctly below Tc, the double-peak distribution t(^) exhibits a deep 
minimum for densities p in between the two coexisting phases pioL [vapor] and pcoL [liquid] .— 
Consequently, a system starting with a low vapor-like density would hardly ever make the 
transition to the liquid-like state or vice versa. Hence, the relative weights of the two phases 
would not be sampled correctly. This difficulty is overcome by biased sampling methods 
that "drive" the system through the coexistence region such as "multicanonical sampling" ,— 
"Wang-Landau-sampling" ^^ or "successive umbrella sampling"— which has been used in this 
work. In the simplest implementation, the algorithm is constrained to sample configurations 
with only two particles n G (0, 1) in the beginning, and (1, 2) ■ ■ • (n-l,n) later on, span- 
ning the relevant range of densities. The probability distribution can then be calculated 
recursively: 

= Hi^oH2,i ■ ■ ■ Hn,n-1 , (13) 

with Hj j_i being the frequency of occurrence of the j^^ particle over the frequency of occur- 
rence of the (j — 1)**^ particle in the sampling of the (j — 1, j) window. For a more detailed 
description of this method and its extension we refer to Virnau et al.—^'^. Biased grand 
canonical methods have the additional advantage that the minimum in P^ xin) at densities 
near the density of the rectilinear diameter pd{T) 

Pd(r) = (piL + pil)/2 (14) 

is also sampled rather accurately. This minimunt^i^i^i^ corresponds to a free energy 
barrier AF ^ 2'j{T)L'^ which arises from the formation of two (planar) vapor-liquid in- 
terfaces of area L^, each connected with itself via periodic boundary conditions. In this 
expression, 7(T) is the vapor-liquid interfacial tension. For p near Pd{T), the system is in a 
state of two-phase coexistence, a slab-like liquid domain is separated from the vapor via those 
interfaces. Coexisting gas and liquid phases have the same free energy. Therefore, AF is the 
free energy of the interface. It has been amply verified for a variety of systems^'^'^i^'^i^ 
that the relation^ 

PM,T(nd)/P^,T(ncoex) oc exp[-27(T)LVfcBT] (15) 
(where = pd{T)L^ and ncocx = PmclL^) is a valid description of the simulation results. 
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and can be used to extract rather accurate estimates for 7(T). 

Close to Tc the estimates for pcoL, PcoL, pd, and 7(T) suffer from systematic finite size 
effects. It turns out, however, that the finite size effects for are numerically rather small. 
Therefore, the critical density pc can be estimated from pc = Pd{Tc) with Eq. f[T^ . pcoL and 
Pcoex are just the peak values of the density resulting from the equal area rule at Tc. (We 
note that pcoex(T'c) > pcoex(Tc) for any finite L. The peak values only merge into a single 
point Pc at Tc in the thermodynamic limit.) 

The behavior of the density near the critical point can then be obtained, too. In the 
critical region the critical exponent /3 has to take the value /3 = 0.325 of the Ising model 
universality class^ 



pioL - P<iiT) = -B{1 - T/T^r 

P^l.-P^iT) = +B{l-T/T^r. (16) 

Here, the critical amplitude B can be estimated by fitting the actual simulation data in 
the range 0.02 < 1 — T/Tc < 0.1 to Eq. (fT6l) . Note that the left boundary of this interval 
is chosen such that for the typical linear dimensions, finite size effects on the peak position 
estimates for pcoL, PcoL are still very small. The right boundary of the interval is chosen in 
order to justify the neglect of correction terms to the leading term written in Eq. (fT6|l which 
only describes the asymptotic behavior in the limit^ 1 — T/Tc — > 0. 

Our data for the coexistence curve and interfacial tension were derived from an elongated 
box L X L X 2L with size L = 9cr and L = 6.74cr (the latter only very far from the critical 
point). The critical points (Figs. [H [2]) were computed using cubic boxes of size 9a and 
11.3 a. In a few larger box L = 13.5 a was implemented to check the finite size 

effects. After coexistence densities were determined, simulations at coexistence gas density 
were carried out in the NVT ensemble to obtain the coexistence pressure from the standard 
virial expression. 
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III. NUMERICAL RESULTS FOR CARBON DIOXIDE: COMPARISON WITH 
EXPERIMENT AND SIMULATIONS OF ATOMISTIC MODELS 



Figs. [MS] present the coexistence curve, the vapor pressure at coexistence, and the inter- 
facial tension as a function of temperature, and compare them to pertinent experimental 
data.^^ If quadrupolar interactions are neglected {qc = 0), a distinct discrepancy between 
the experimental data and the simulations can be observed for the liquid branch of the 
coexistence curve.-i^ Agreement with experiments improves considerably for the isotropic 
quadrupolar model. A value of qc = 0.387 was used which corresponds to the experimental 
value of the CO2 quadrupolar moment \Q\ = 4.3 DA (Eq. flTU]) ) as discussed above. It is also 
very gratifying that both coexistence pressure (Fig. H]) and interfacial tension (Fig. are in 
almost perfect quantitative agreement with experimental data, although for these quantities 
there is no adjustable parameter available whatsoever. In particular, the interfacial tension 
for qc = deviates from the experimental data rather distinctly, while for = 0.387 there 
is excellent agreement. 

A small but systematic discrepancy is still present for the liquid branch of the coexistence 
curve (Fig. [3]). Hence, we have also tried to take qc as an adjustable parameter to optimize 
the agreement between the simulated coexistence curve and the experiments. The rationale 
for doing so is twofold: first, there is a considerable uncertainty in the experimental value 
for Q, leading to a 30% uncertainty in - it is not even clear that the value of Q for CO2 
in the vapor phase and in the liquid are exactly the same. Secondly, it might be better to 
choose an effective value for Q because our spherically symmetric model (Eq. ([5])) is a rather 
incomplete description for the interactions between elongated CO2 molecules. In principle, 
the systematic coarse-graining of a chemically realistic model could lead to some effective 
value for Q, which is larger than the experimental one. 

Thus, Figs. [3ll5] also include some simulation results for a second choice of gc, namely 
qc = 0.470. Fig. [3] shows that now the agreement between simulation and experiment for the 
liquid branch of the coexistence curve is better than for qc = 0.387, but for the vapor branch 
it is slightly worse. The same slight deterioration of the agreement can also be observed for 
the coexistence pressure (Fig. H]) and the interface tension (Fig. [5]). We conclude that an 
absolutely perfect agreement between any simplified model, such as Eq. ([5]), and a real system 
simply cannot be expected. Some uncertainty about the optimum choice of the parameters 
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of such a coarse-grained model is simply inevitable. Actually, the level of agreement between 
experiment and our model is very good for both choices of qc- This is gratifying, since the 
model will serve as an excellent starting point for the coarse-grained modeling of various 
polymer solutions containing CO2 as a solvent. 

A model of the type of Eq. ^ (named isotropic multipolar or IMP) was also used in 



Ref. 



40 



42 



and the vapor-liquid coexistence curve of CO2 was determined with temperature 
quench MD techniques.— The simulation results of Ref. |421 are reported in Fig. [3] (see o), 
too. Although large systems were used, error bars in the determination of the coexisting 
densities using NEMD are large in comparison with ours as discussed above. (Errors for our 
simulations are smaller than the size of the symbols and therefore not shown in Figs. [3]- El) 



We also note that Ref. 



42 



uses Lennard- Jones parameters that differ significantly from ours, 
namely e/kB = 215.0 K and cr=3.748 A while we use e/kB = 252.8 and cr=3.785 A for 
\Q\ = 4.3 DA. This is mainly related to the larger cutoff radius of 4 cr used in Ref. Q, 
which increases the critical temperature. Our agreement with experimental results (i.e., 
coexistence curve Fig. [3l coexistence pressure Fig. H] and isobar Fig. [8]) is, however, clearly 
very good because our grandcanonical simulations allows for a very precise determination 
of the critical point. 

Let us ask how our simulation results for the coarse-grained model compare to the results 
obtained for atomistic models of CO2. Figs. [6l [7] and [8] present such comparisons for the 
coexistence densities and pressures with some results available in literature. The EPM 
model^ (denoted by + in Figs. [6] and [7]) overestimates the vapor density at coexistence 
and underestimates the coexistence pressure systematically, while the liquid densities are 
underestimated only for T < 260K. For T > 280K, the liquid densities of the atomistic 
simulation are too large due to the overestimation of Tc. When the atomistic model is 
rescaled (EPM2)^ so that the critical temperature and density are matched (denoted by < 
in Figs. [6]and[7j), the agreement between the model calculation and experiment is almost 
as good as for our coarse-grained model. However, the rescaled data for the coexistence 
pressure are slightly but systematically too large. The coexistence line for the EPM2 model 



28l . in agreement with the previous work.^^ In Fig.[8]we include 



has also been obtained in Ref. 
simulation results of Ref. US 



for the EPM2 model for the supercritical isobar (200 bar). The 
both models work very good in the supercritical region, although the coarse grained model 
gives slightly better agreement with experimental data for both choices of qc used in this 
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work. Recently,— another optimized version of the EPM2 model has been proposed in which 
the atomistic energies, lengths and charges have been rescaled to optimize agreement with 
the coexistence experiments. As a consequence, the agreement with experimental results 
is very good, in particular for the coexistence pressure (see O in Fig. [7]). Simulations fit 
the experimental curve perfectly below 270K, while for higher temperature small deviations 



appear. In Ref. |33!, two center Lennard- Jones models which include a quadrupolar point 
have been studied extensively, and coexistence densities and pressure were obtained. Tuning 
atomistic parameters, the agreement with the experimental curve has been optimized^i 
without any physical input. As a result, a quadrupolar moment for CO2 predicted in Ref. 



34 



equals |Q|=3.7938 DA which is quite off from the experimental value 4.3 DA. Finally, 
there is also a recent simulation, which uses two ab-initio potentials named BBV^^ (denoted 
by □ in Figs. [6] and [7]) and SAPT-s^^ (denoted by o in Figs. [6] and [7]). Results are quite 
off the respective experimental values, but unlike to the previously mentioned models, no 
fitting procedures have been applied. No data on the interfacial free energy of the atomistic 
model are available so far to which we could compare our results. Figs. [6] and [7] demonstrate 
that the rescaled atomistic model agrees better with experiment than the simple LJ model 
which ignores the quadrupolar interaction completely. However, in comparison with the 
present model (Eq. ([5])), the atomistic models offer no advantages, even if one rescales the 
parameters to match the critical point. In fact, the use of Coulomb interactions in the 
atomistic models makes the code considerably slower. 



IV. OTHER QUADRUPOLAR FLUIDS 

For a detailed discussion on how to derive simulation parameters for an arbitrary 
quadrupolar substance, the reader is referred to appendix [Bl Here we would like to fo- 
cus on testing the model for other quadrupolar substances. Using literature data for Q, 
Tc^exp and Pc,exp for various molecular fluids, we can use our master curves (Fig. [1]) to predict 
the value of qc and describe these fluids with our model, Eq. ([5]). Inserting Eq. ([8]) into Eq. 
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we obtain 



-1 10/3 



Qc = 



PcjCxp 




= X. 




(17) 



Note that Aexp contains all the experimental parameters which are required to define the 
model. Fig. [9] plots qc as a function of Acxp for CS2, N2, CO2, C2H2, and CeHg. 

One recognizes immediately that for N2 and CS2 the effects of the quadrupolar interac- 
tions can only be minor, since qc is very small. Consequently, the simple LJ model (where 
quadrupolar effects are completely neglected) should be a reasonable description of the co- 
existence densities, coexistence pressures, and interfacial free energies of those fluids. Fixing 
the LJ parameters for N2 via Tc and pc as done in our previous work,-i^ we can test im- 
mediately this hypothesis (Fig. [10]). As expected, the deviations from the simple LJ fluid 
are indeed much less pronounced than for CO2. Note that these deviations between the 
measured and the predicted coexistence curves for these fluids with small qc are comparable 
to the deviations found between the simple Lennard- Jones coexistence curve and the exper- 
imental results for noble gases such as Ne, Ar, Kr and Xe. These systems are considered to 
be the best experimental realization of a Lennard- Jones fluid (Fig. [TTl) . In a rescaled rep- 
resentation (T/Tc plotted vs. p/pc), however, the various noble gases do not exactly satisfy 
a "law of corresponding states" . This implies that even for systems with perfectly spherical 
atoms, a description in terms of (classical) point particles interacting with purely pairwise 
potentials of the same functional form (with one parameter for the strength and another for 
the range) is not strictly valid. 

These small deviations may be due to the need for three-body forces^, or quantum cor- 
rections which account for differences in atomic masses. The inclusion of the three body 
interaction is computationally extremely expensive. Indeed in the evaluation of the total en- 
ergy of the system one would need to evaluate a total number of contributions that scales like 

instead of N^ as for the two body interactions (N, being the total number of molecules). 
For this reason the inclusion of such effects in our simple (and cheap) modeling is out of 
discussion, especially in view of more complicated polymer solution applications. There are 
several attempts^"^ which try to capture the three body interaction in an effective (density 
dependent) two body interaction. These methods cannot be used in non homogeneous fluids 



16 



and generally where strong density fluctuations are present, like near the critical point. The 
fact that the method proposed in this work is based on a careful investigation of the critical 
points of the coarse grained models invalidates the scheme proposed in Ref. 
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69 



However 



in Ref. 



671 a quantitative estimate of the effects of the three body interaction is given starting 



from a careful scaling investigation of the rectilinear diameter (11 4p 

with a ~ 0.11. The authors shows that in Eq. flTSl) A^^^ is related to the fleld mixing effect 
(indeed the lack of the particle hole symmetry), while could give an estimate of the three 
body interaction. A Mean Field van der Waals equation predicts^ A\ = 2/5. Deviations of 
the experimental data from this law of corresponding states {Ai = 2/5) are supposed to be 
related to the emergence of another energy scale like that of three body interactions. Fig. 4 



of Ref. 



671 suggests (for CO2) Ai ^ 0.95 which differs signiflcantly from the van der Waals 



value Ai = 0.4 but is comparable with other fluids in particular Xenon. Comparing now the 
predictions for Xenon (Fig. [TT]) and Carbon Dioxide (Fig. [3]), one can easily conclude that in 
our case the quadrupolar interactions are much more relevant than three body interactions. 

For the sake of completeness, in Fig. [11] we have also included the full LJ potential. In an 
unsealed representation one would of course observe large differences between the results for 
the full Lennard- Jones potential and those for the cut-and-shift Lennard- Jones potential. In 
a scaled representation these differences vanish almost completely (except for small densities 
on the gas branch of the binodal) so that due to its computational efficiency the cut-and-shift 
potential should be preferred in coarse-grained simulations. 

The case of benzene (CeHg) is even more interesting. Depending on which experimental 
value is adopted for Q, one flnds Qc in the range from = 0.121 (for Q = 10 DA) to 
Qc = 0.247 (for Q = 12 DA). Fig. [T2] compares experimental values for the coexistence 
densities, coexistence pressure and interfacial tension with our predictions, using Qc = 0.247. 
In this case we also observe a clear improvement of the agreement with experimental data 
with respect to the pure Lennard Jones case {qc = in Fig. [T2l) . Deviations are of the same 
order of magnitude as for nitrogen (Fig. [101) and noble gases (Fig. [TT]). 
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V. PREDICTIONS FOR THE EQUATION OF STATE RESULTING FROM PER- 
TURBATION THEORY (PT) 



In this section we present results for coexistence densities and coexistence pressures 
(Fig. [T^ . which were obtained analytically using an equation of state™ in the Mean Spher- 
ical Approximation (PT-MSA)™. As is well-known^^ such approaches should work well at 
temperatures and densities away from the critical region.^-^ This expectation is reconfirmed 
by our results (Fig. [T3l) . which show good agreement at temperatures below 0.9 T^. For 
0.9 Tc < T < 1.2 Tc, there are distinct deviations between simulations and theory because 
PT-MSA overestimates the critical temperature by about 10% and furthermore the slope 
of the binodal in the critical region is mean-field-like in PT-MSA and Ising-like in the sim- 
ulation. For low temperatures the deviations are quantitatively smaller, however, the MC 
results and PT-MSA results cross at T* ^ 1 (if = 0.387) and T* ^ 1.05 (if = 0.470) on 
the liquid branch. Note that our comparison involves no adjustable parameter whatsoever. 
For many practical applications one will be interested in the temperatures and/or densities 
outside the critical region. Hence, the results shown in Fig. [13] are encouraging in that a 
relatively simple analytic method such as PT-MSA (see Appendix A for some details on this 
method) works well as a description of the equation of state for molecular fluids like CO2 
away from the critical region if an isotropic quadrupolar interaction is included. To some 
extent this minimizes the need for massive Monte Carlo (MC) efforts to explore phase space. 
Even though MC simulations are required to determine e, a and qc from Tc^exp, Pc,exp and 
Q, the results are already contained in Fig. [Hand Eqs. ( I9all9bl) . Therefore, no new efforts 
with MC simulations will be needed for any future applications of PT-MSA in the context 
of our model. 



VI. CONCLUSIONS 



In the present work, the thermodynamic properties of a coarse-grained model for 
quadrupolar fluids were investigated. A particular emphasis was put on the question to 
which extent the equation of state and the interfacial tension between coexisting vapor and 
liquid phases can be described accurately. 

The aim of this work hence is not a chemically detailed modeling of quadrupolar fluids 
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on an atomistic level, but rather to derive a model which is bot simple and accurate enough 
that it can serve as a starting point for the description of binary fluid mixture, solvents 
in polymer solutions, etc.. Obtaining efficient models for such purposes is a topic of great 
current interest. 

As experimental input parameters, our description only requires knowledge of the ex- 
perimental critical temperature Tc^exp and the critical density pc,cxp of the fluid and the ex- 
perimental quadrupole moment Q of the molecule. The quadrupolar interaction is treated 
in a spherical approximation^!^ which can be derived from thermodynamic perturbation 
theory. This leads to an effective potential proportional to Q'^/{Trlj), where T denotes the 
temperature and r^j the distance between the centers of mass of molecules i and j. The 
application of the isotropic quadrupolar interaction is mainly motivated by the desire to 
have a very fast simulation code. Steric and dispersion forces are simply modeled by a 
Lennard- Jones potential involving parameters e and a, which define the strength and the 
range of the interaction, respectively. In practice, the potential is cut and shifted to zero at 
a cutoff range Tc = 2v^, which is again motivated by our desire to speed up calculations. 
We also provide evidence that this particular approximation mostly affects the conversion 
factor from e to experimental temperature and hence does not alter results significantly. 

For the description of a real system, simulation parameters e, a and qc = 
Q'^/{scr^'^kBTc^cxp) need to be determined from experimental values Tc^exp, Pc,exp and Q in 
physical units. To address this problem, we have determined master curves T*{qc) /T*{0) 
and P*{Qc)/Pc{^) as a function of qc (Fig. [1], Eqs. ( I9al) .( l9bl) ). This task is performed eas- 
ily using grand-canonical Monte Carlo simulations^'^'^ in combination with reweighting, 
successive umbrella sampling^ and finite-size scaling methods^"^"^. With modest computa- 
tional effort, these master curves are determined with a relative accuracy which is distinctly 
better than 1%. 

Carbon dioxide is a prototype of a linear elongated molecule with a rather large 
quadrupole moment. Comparing our predictions for the coexistence curve, vapor pressure 
at coexistence and interfacial tension with corresponding experimental data^^, we found 
encouragingly good agreement (Figs. [3lHf5|) . Note that after having fixed the scales for tem- 
perature and density via e and a, no further parameters need to be adjusted, neither for the 
pressure (FigH]), nor for the interfacial tension (Fig. [5]). The level of agreement which we 
have achieved is clearly nontrivial. However, the inclusion of quadrupolar effects is essential 
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to the model and agreement with experiments deteriorates significantly if CO2 is described 
by a Lennard- Jones particle without quadrupole moment. 

Our model produces rather accurate off-critical isotherms, too. As expected, the compar- 
isons also reveal small discrepancies, since such a simple model cannot be absolutely perfect. 
However, a more realistic model, based on an all atom description of CO2 which involves 
considerably more complicated potentials, performs distinctly worse in comparison to our 
model - except if experimental critical parameters are used to empirically re-calibrate the 
atomistic potential. In our view, such a procedure looses the advantage of a fully predictive 
modehng that does not need experimental input. Complicated atomistic models also lead 
to rather slow simulation programs (partial charges require to deal with rather long range 
coulombic interactions, etc.). While such models may still be manageable for the simulation 
of pure fluids, their drawbacks become clearly apparent when the approach is extended to 
binary or ternary fluids. In mixtures, a large control parameter space needs to be explored 
and several phase separations may compete with each other, leading to very involved phase 
diagrams. 

We emphasize that our successful description of carbon dioxide is by no means accidental. 
As a counterpart, we also consider nitrogen, a fluid with a considerably smaller quadrupole 
moment. In this case, a simple Lennard- Jones model with no quadrupolar forces should 
provide an equally good description, and in fact it does. The deviations are comparable 
to the deviations found between the coexistence curve of " Lennard- Jonesium" and those of 
various noble gases (that do not superimpose precisely in a re-scaled representation shown 
in Fig. [H] either.) This indicates that a simple pair potential with two parameters for the 
scales of energy and range does not suffice even for these prototypes of simple spherical 
atoms. 

As a further example, we also present a comparison between our model and experimental 
data for benzene (CeHg). Again, the agreement is very good. This result is of great interest, 
since the shape of the benzene molecule differs considerably from CO2, consisting of a disk 
rather than an elongated ellipsoid. 

A very interesting question is the extent to which this concept can actually be carried 
over from simple fluids to binary mixtures and polymer- solvent systems. Are interactions 
between different types of molecules captured by simple Lorentz-Berthelot mixing rules, 
when one describes the pure constituents with the quality of the present work? We shall 
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address this very interesting and potentially practically useful question in a forthcoming 
paper. We also hope that the present work will stimulate some analytical research, starting 
from general statistical mechanics of fluids, to provide a better theoretical understanding for 
the high accuracy of our approach. We also point out that the knowledge of the appropriate 
parameters e, a and qc allows a rather accurate description of the equation of state by liquid- 
state perturbation theories at state points sufficiently away from the critical region (Sec. V). 
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APPENDIX A: MEAN SPHERICAL APPROXIMATION (MSA) PREDICTIONS 

In this appendix we want to give some technical details concerning the analytical pre- 
dictions presented in this paper. For more details we refer to the original literature. In 
particular, the equation of state (EOS) used in this work is a straightforward generalization 



of the EOS given in appendix B of Ref. 



70 for the case in which four Yukawa tails are used 
in which the Ornstein-Zernike (OZ) 



72 



73 



instead of two. We follow the strategies of Refs. 
equation is solved in a first order MSA closure. The general ideaA is to divide the potential 
into a repulsive part (that becomes the reference potential) plus a perturbative attractive 
part 

Ux{r) = ^ (Al) 
[ >^Uatt{r) if (To < r < rent, 

where U{(To) = 0, U^epir) > 0, f/att('") < and A is the perturbative parameter. The 
reference system (A = 0) is modeled by hard spheres with a proper radius d^s,— computed 
using t/rep— • In order to get corrections to the reference free energy Aj-^f, a systematic 
expansion in A is developed (the general expression for A — A^-ei is standard and can be 
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found for example in Ref. 



been obtained in. Refs. 



72 



70 
73 



(Eq. B5)). The explicit solution up to second order in A has 



The key point developed in Ref. 



73l is to fit Uatt with a couple 



of Yukawa tails. In the case of the LJ potential this yields 



ttLJ ^ 



g-2;i(r-(To) Q-Z2{r-cro) 
-Ci h C2- 



r 



r 



(A2) 



In this work the LJ part of the potential is fitted using the same Yukawa tail as reported 
in Ref. |70| (Eq. B6). Equation (1A2I) allows us to invert some Laplace transforms that are 
present in the Tang-Lu solution''^ and to obtain an analytical expression for the free energy 
which is explicitly given in Eq. B7-B10 of Ref. [70 for the apolar-fiuid case g = 0. 

For the general case q ^ 0, Eq. ( lAll) will induce the same decomposition on both the LJ 
part and quadrupolar part of the potential 



att(rep) 



u. 



LJ 
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att{rcp) a'tt{i'cp)' 



(A3) 



In ( ]A3I) we have used two more Yukawa tails to fit the quadrupolar interaction Ul^^ 



U. 



IQQ 



g-23(T--CTo) g-24(r-(To) 

-C3 h C4- 



r 



r 



3^QQ(Q,z„ao;r). 



(A4) 



Because q (= qcTc/T) is factored out in flA3p . 03^4 and 2:3^4 do not depend on temperature T. 
This is an important simplification because using (]A2p and (]A4I) we can get an immediate 
fit for f/att f!A3p for every q and T 
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att 
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(A5) 

to the case in which 



By using the previous fit (lASP and extending Eq. B7-B10 in Ref. 
more than two Yukawa expressions are used to fit the potential, we have obtained the desired 
EOS used in the present work. 



APPENDIX B: DETERMINATION OF SIMULATION PARAMETERS 

Simulation parameters e, a and qc are needed to convert simulation units into experi- 
mental units. Knowledge of qc, or rather q = qc ■ Tc/T is also required as input before a 
simulation can be started. In Table [H we have collected the simulation parameters for the 
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quadrupolar substances mentioned in the paper. However, we would also like to convey 
some hands-on knowledge on how to calculate these parameters and extend the model to 
substances not listed in Table [H Furthermore, we provide fitting curves (Table [ITl) which 
allow us to determine the phase diagram of an arbitrary substance without additional MC 
simulations. 

For qc=0, e and a can be determined directly from the critical temperature Tc and the 
critical density pc using Eq. ([8]). For qc ^ 0, the location of the critical point itself depends 
on qc. Therefore, e and a also depend on (Eq- dHD), and a simple iteration procedure 
can be formulated. Starting with gc=0, Tc and pc are computed using the master curves 
from Eqs. fl9al) and (]9bp . From these results, e and a are determined with Eq. ([H]) and a 
new value for qc with Eq. The iteration is repeated until q^ e and a converge. Usually, 
around 5-10 iterations are sufficient to obtain simulation parameters with good accuracy 
without any additional simulations. In the following, we present a pseudo-code for our 
CO2 calculations which can be extended to any quadrupolar substance by substituting 
experimental values for Q=4.3 DA, Te=304.1282 K, and pc=10.6249 mol/l: 



Initialize variables 

Q = 4.3 /* Dk */ 

Q = Q*3.33564*10-^°; /* convert Q to SI units */ 

Tc,exp = 304.1282; /* K */ 

rhoccxp = 10.6249; /* mol/l */ 

TLj(q=0) = 0.99821 /* critical temperature of simulation for q=0 */ 

rhoLj(q=0) = 0.32276 /* critical density of simulation for q=0 */ 
q = 0; 



Iteration 

for (i=0;i<20;i++) { 

T = TLj(q=0) * (1 + 0.46111 * q + 0.17571 * q^); /* Eq.(l9aD */ 

density = rhoLj(q=0) * (1 + 0.19298 * q); /* Eq.(l9b]) */ 

epsilon = Tc,exp * 1.38065 * lO^^s / T; /* Eq.® */ 

Sigma = (rhocexp * 1000 * 6.02214 * lO^^ / density ) "^/^ ; /* Eq.® */ 

Qi = Q/(sqrt(epsilon*sigma^)); /* Eq.® */ 
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q = Q4 / (T * 1.237990147 * IQ-^O); /* Tsi„,=kBTe.p/£, (47r£o)' - SI units */ 
print T, epsilon, sigma, q; 

} 

Alternatively, qc can also be determined from the fitting curve in Fig. [TDl Agxp is a 
dimensionless parameter, which already contains all the experimental information required 
to define the model. If all constants are included, Agxp reduces to 

Aexp = 96.754 ■ 10-5-^(p,,,xp)f . (Bl) 

c,oxp 

In this equation, one simply needs to plug in experimental values for quadrupolar moment 
Q in DA, critical temperature T^^oxp in K, and critical molar density Pc,exp in mol/cm'^. q^ 
can be read off from Fig. [10] or determined via the following fit to the curve: 



qc = Aexp(43.1018 - 266.251Aexp + 5047.01A^^p) Aexp < 0.02 (B2) 

Tc, Pc £ and cr follow from Eqs. fl9aj) . fl9b|) . and ([8]). 

Finally, we demonstrate how our accumulated simulation data can be used to provide a 
rough estimate of the phase diagram for an arbitrary quadrupolar substance without any 
additional MC simulations. We simulated several values for qc in the range of 0.1 < gc < 0.47. 
Four temperatures were considered such that T*{qc)/T*{qc) (i = 1, ■ ■ ■ , 4) is independent of 
q^. T* = 0.974499 ■ T*, T* = 0.932125 ■ T*, T* = 0.864337 ■ T*, and T; = 0.813494 ■ T*. 
As indicated before, critical quantities scale almost linearly with qc (Fig^H Eqs. (pal) and 
(I9b|) ). During our investigations, we observed that this approximation also holds away from 
criticality. The corresponding fitting curves are listed in Table [TTl 

First, one needs to determine e, a and qc for the substance in question as demonstrated in 
the previous section. Vapor and liquid coexistence densities, interface tension and pressure 
at the selected temperature can be computed by inserting qc into the respective fitting 
curves. The following equations can be used to convert the results from simulation units to 
experimental units: 

_ g(gc) _ * _ , g(gc) _ . g(gc) 
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TABLE I: Experimental data and simulation parameters for several quadrupolar substances as 
obtained in the present work 



TABLE IL Fitting curves to determine coexistence properties for an arbitrary quadrupolar sub- 
stance at selected temperatures T* (qc) /T* [qc) (i = 1, • • • , 4): = 0.974499-Tc*, = 0.932125-r*, 
Tg* = 0.864337 • T* , and = 0.813494 • T* (see text) 

FIG. 1: Master curves: normalized critical temperature T*{qc)/T*{0), normalized critical density 
Pcilc)/ pI{0), and normalized critical pressure PciQc) / Pci^) plotted versus the quadrupolar param- 
eter Qc- Symbols represent simulation data, curves are the interpolating functions (Eqs. ([9a]) and 
&h\i ) and p*(gc)/Pc(0) = (1 + 0.67423^^ + 0.274349 g^) with p*(0) = 0.08722L 



FIG. 2: Second and fourth order cumulants U2, U4 plotted for q = 0.3 versus T* = ksT/e for 
three choices of L. Broken horizontal values indicate the theoretical values established for the 
Ising universality class.— 1^ From the intersections one can conclude T* = 1.152 it 0.003 for this 
particular case. Inset: the slope of the fourth order cumulants (Yi) as a function of the box size, 
on a log-log scale. The data points fall on a straight line with a slope equal to 1.584 in agreement 
with the finite size prediction l/u, with u ~ 0.630 for the Ising universality classj^ 

FIG. 3: Coexistence curve of CO2 plotted in the temperature-density plane. The broken curve 
denotes the experimental data (from NIST— ), the full curve is the result for the LJ model without 
quadrupolar interactions^. Solid square denotes the critical point of CO2. (x) and (*) are the 
results of the present pVT work for two choices of Qc = q(Tr) as indicated in the figure, (o) are the 



results of the spherical averaged model investigated in Ref. 



42. 



FIG. 4: Coexistence pressure of CO2 plotted vs. temperature. The broken curve denotes the 
experimental data,~ the full curve: the results for the LJ model without quadrupolar interactions, 
(x) and (*) are the results of the present NVT work for two choices of qc = q{Tc) as indicated in 
the figure. 
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FIG. 5: Interface tension ^{T) of CO2 plotted vs. temperature. The broken curve denotes the 
experimental data (from NISTSS), the full curve: the results for the LJ model without quadrupolar 
interactions.^ (x) and (*) are the results of the present work for two choices of Qc = q{Tc) as 
indicated in the figure. 



FIG. 6: Coexistence curve of CO2 plotted in the temperature-density plane. The broken curve 
denotes the experimental data (from NIST— ), the full curve: the results for LJ model without 
quadrupolar interactions^. (•) denotes the critical point of CO2. (*) and (x) denote the results 
of this work for Qc = 0.387 and qc = 0.470, respectively. (+) are the results of the EPM model 



introduced in Ref. 



271 . (v) are the results from Ref. 



271 for the EPM model with flexible molecules, 



which give essentially the same thermodynamic properties as the rigid molecules. (<) are the results 



of Ref. I27I for the rescaled EPM model (EPM2). (o) and (□) correspond to simulations^^ of two 
ab initio potentials.—*^ 



FIG. 7: Coexistence pressure of CO2 plotted vs. temperature. Labeling of curves and symbols is 
the same as in Fig. O We also show simulations of an optimized EPM2 model^S (see O) which 
is in good agreement with experiments. We stress that the nice agreement of our model with 
experiments near the critical point is not given a priory because our method only fixes the critical 
temperature and the critical density. 



FIG. 8: Supercritical isobar for p=200 bar. The broken curve denotes the experimental data^^ 
(x) and (*) are the results of the present NVT work for two choices of Qc = |(7 );) as indicated in 
the figure. (<l) are the prediction of the atomistic EPM2 model given in Ref. 2^. The coexistence 
curve near the critical point is also reported. 



FIG. 9: Estimates for the quadrupolar parameter qc for various quadrupolar fluids characterized 
by parameter Agxp (Eq. (jl7p ). The corresponding experimentally measured quadrupole moments 
Q of these systems are quoted in brackets (see also Table H]) . 
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FIG. 10: Lennard- Jones results {qc = 0) for N2. From top to bottom: coexistence cm've in the 
temperatm'e-density plane, vapor pressure vs. temperature and interface tension vs. temperature. 
Symbols correspond to simulations of a simple Lennard-Jones model without quadrupolar moment 
obtained from /iVT simulations^ (coexistence densities and interface tensions) and NVT simulations 
(pressure). The broken curves denote the experimental data (from NIST— ). 



FIG. 11: Coexistence curves (T/Tc plotted vs. p/pc) for various noble gases in comparison with the 
prediction of the cut-and-shifted Lennard-Jones model (LJ)^ and the full Lennard-Jones model.— 



FIG. 12: New predictions for benzene (CgHe). From top to bottom: coexistence curve in the 
temperature-density plane, vapor pressure vs. temperature and interface tension vs. temperature. 
The broken curves denote the experimental data,— the full curve is the result of the simple Lennard- 
Jones model. (>) denote the present results which include an isotropic quadrupolar interaction for 
Qc = q{Tc) corresponding to Q = 12 DA. 



FIG. 13: Coexistence densities and coexistence vapor pressure: a comparison between the MC 
simulations and the PT-MSA prediction. The two choices of qc = q{Tc) used in this work are 
included as indicated. 
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Subst. 


Q [dA] 


-^c,cxp [K] 


Pc,exp [mol/1] 


'^cxp 


Qc 


e/kB [K] 




CO2 


4.3 


304.1282 


10.6249 


0.009430 


0.387 


252.829 


3.785 


CS2 


3.6 


552 


5.78 


0.0001848 


0.0080 


550.95 


4.528 


N2 


1.47 



126.2 
126.2 


11.18 
11.18 


0.0008864 



0.038 



124.208 
126.426 


3.642 
3.633 


C2H2 


5.5 


308.3 


8.913 


0.013775 


0.553 


235.942 


4.052 


CqB.q 


12 


562 


3.9 


0.0059311 


0.247 


500.468 


5.242 



TABLE I 
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Observable 


Fitting Formula 


* 


0.162 


Pin 


0.099506 -0.0094827^0 


' *-*iy 


0.055372 -0.017106^0 


P\a 


0.036003- 0.018 


* 


0.49215 + 0.12426 qc + 0.021146 ql 


Ph 


0.57055 + 0.15313 q^ + 0.025081 ql 


* 


0.64597 + 0.1531 qc + 0.09854 ql 




0.68355 + 0.22094 gc + 0.042765 ql 


/ i 


0.020384 + 0.016672 q^ + 0.027991 ql 


72 


0.068376 + 0.072064 + 0.082864 


73 


0.16187 + 0.19493 qc + 0.18704 ql 


74 


0.23945 + 0.29931 qc + 0.30352 ql 


P^ 

^ 1. 


0.075861 + 0.041526 + 0.024072 


P*2 


0.056804 + 0.026873 + 0.011408 g^ 


Pi 


0.035115 + 0.0099939 gc + 0.00067637 


Pi 


0.023617 + 0.0010425 qc - 0.0009939 ql 



TABLE II 
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